The Proposed Quadruple System SZ Herculis: Revised LITE Model and 

Orbital Stability Study 



Tobias Cornelius Hinse 
Korea Astronomy and Space Science Institute, Daejeon 305-348, Republic of Korea 
Armagh Observatory, College Hill, BT61 9DG, Armagh, NI, UK 
tchinseOgmail . com 
Krzysztof Gozdziewski 
Nicolaus Copernicus University, Torun Centre for Astronomy, PL-87-100 Torun, Poland 

Jae Woo Lee 

. Korea Astronomy and Space Science Institute, Daejeon 305-348, Republic of Korea 

Nader Haghighipour 

Institute for Astronomy & NASA Astrobiology Institute, University of Hawaii, 96822 HI, USA. 

Chung-Uk Lee 

Korea Astronomy and Space Science Institute, Daejeon 305-348, Republic of Korea 

ABSTRACT 

In a recent study, Lee et al. presented new photometric follow-up timing observa- 
tions of the semi-detached binary system SZ Herculis and proposed the existence of 
two hierarchical cirumbinary companions. Based on the light-travel time effect, the 
two low-mass M-dwarf companions are found to orbit the binary pair on moderate to 
high eccentric orbits. The derived periods of these two companions are close to a 2:1 
mean-motion orbital resonance. We have studied the stability of the system using the 
osculating orbital elements as presented by Lee et al. Results indicate an orbit-crossing 
architecture exhibiting short-term dynamical instabilities leading to the escape of one of 
the proposed companions. We have examined the system's underlying model parameter- 
space by following a Monte Carlo approach and found an improved fit to the timing 
data. A study of the stability of our best-fitting orbits also indicates that the proposed 
system is generally unstable. If the observed anomalous timing variations of the binary 
period is due to additional circumbinary companions, then the resulting system should 
exhibit a long-term stable orbital configuration much different from the orbits suggested 
by Lee et al. We, therefore, suggest that based on Newtonian-dynamical considerations. 
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the proposed quadruple system cannot exist. To uncover the true nature of the observed 
period variations of this system, we recommend future photometric follow-up observa- 
tions that could further constrain eclipse-timing variations and/or refine light-travel 
time models. 

Subject headings: binaries: close — binaries: eclipsing — stars: individual (SZ Herculis) 



1. INTRODUCTION 

Formation of multipl e star systems is complex and is believed to occur either by interac- 
tion/capture mechanisms (jvan den Berk et al.l 120071 . and references therein) during the formation 
and dynamical evolution of globular star clusters, or directly from , a m a ssive primordia l disk 
involving accretion processes or disk instabilities (|Lim &: Takakuwal 120061 : iDuchene et al.l 120071 : 
Marzari et al.ll2009l ). The exact formation channel is not yet fully understood. However, these 
mechanisms might be capable of producing v arious star systems characterized by different orbital 
architectures and/or hierarchies (|Evanslll968l ). 



Cirumbinary objects belong to a category of hierarchical star systems where one (or more) 
massive companion(s) orbits a pair of stars. An example of such systems is the quadruple (or 
quaternary) system of HP 98800 consis ting of two distinct spectroscopic binaries orbiting a com- 
mon center of mass (iFurlan et al.l 120071 ). Single or multiple low-mass circumbinar y companions of 



planetary nature have been recently discovered frorn grou nd-based observations ()Lee et al. 
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(NN Ser, HW Vi r, HU Aqr and DP L eo) have proposed orbital properties that seem to render their 



orbits unstable ( 
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In a recent work, iLee et al.l (j201ll ) presented new photometric observations of the Algol-type 
semi-detached binary star system SZ Herculis (SZ Her ( AB) hereafter). Based on more than 1,000 
eclipse measurements (spanning more than a century) these authors were able to detect a significant 
change in the systems orbital period manifesting itself as eclipse timing variations (ETVs). Such 
a change in the binary period can be due to i) the interaction between the two binary components 
through their magnetic fields, mass-transfer, or tidal interactions (resulting in apsidal motion) ii) 
gravitational perturbations by additional massive obj ects (companions) and/or iii) by the light- 
travel time effect (LITeI^) also known as R0mer effect (Irwin 1952 . 19591 ). 



It is important to note that timing measurement errors can be uncorrelated (white noise 



^this notation follows the notation as suggested by Hessman et al. ( 201oh 
^Sometimes referred to as LTTE or LTT in the literature 
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following a Gaussian distribution) and/or of systematical (correlated o r red noise) origin . For 
timing measurements of pulsars, this has recently been pointed out by I Coles et al.l (|201ll ) as a 
possible cause of errors in estimating the model parameters. In the past, neglecting the eff ect of 
red noise was responsible for the false detection of planets around pulsars (jBailes et al.lll99ll ). 



The LITE effect implies the presence of one or more massive object (s) which can result in the 
reflex motion of the binary barycenter about the total system's center of mass. This reflex motion 
(or binary wobble) gives rise to time varying light-travel paths resulting in differences in the periodic 
mid-eclipse timing. It is important to note that the LITE effect is a geometrical effect and does 
not involve gravitational perturbations. For instance, in a hierarchical triple stellar system (a third 
companion orbiting a close binary), the third object crea tes a single -body LITE effect, introducing 
a sinusoidal-like variation in the binaries orbital period (jlrwinlll959l ). 



In their work, iLee et al.l (j201ll ) applied various eclipse timing variation models in an attempt 
to describe the observed period variations of SZ Her. Their most promising fit to the times of 
minimum light points to a two-companion LITE model, implying the existence of two low-mass 
M-type stars orbiting SZ Her on circumbinary orbits. The orbital periods of these companions are 
P3 86 and P4 43 years, suggesting a near 2:1 mean- motion resonance between their orbits. 



In this work, using the fitted orbital parameters from lLee et al.l ( 201ll ). we study the stability 



of the two suggested M-dwarf companions. We then perform an extensive parameter-search for a 
best-fit model using the complete set of timing data of SZ Her as compiled from the literature and 
transformed to the terrestrial time (TT) standard. In particular, we carry out a quasi-global Monte 
Carlo search of a variety of two-companion LITE models to explore a region of parameter-space. 
Finally, the orbital stability of our best-fit model is examined and its implications are discussed. 



2. BARYCENTRIC TWO-COMPANION LITE MODEL 

The effect of a single circumbinary companion on the binary period has been presented in llrwi: 



(|l952l ). In that study, the eclipsing binary pair is regarded as a single object with a mass equal 
to the sum of the masses of each component. Because of the presence of an additional massive 
body (e.g. a stellar or planetary companion) the single binary object follows an orbit around the 
system's barycenter. This orbit, known as the light-travel time orbit or LITE orbit, gives rise to a 
modulation of the eclipsing period due to changes in the light-travel distance. We r eproduce the 



LITE orbit (solid ellipse) in Fig. [T] along with the orbit of the additional companion. Ilrwini (jl952l ) 
presents a discussion of the light-travel time orbit in a coordinate system with origin at the center 
of the LITE orbit. For reasons of symmetry, when discussing properties of the LITE orbit, this 
choice is suitable for the geometric interpretation of the resulting light-travel time curve (or O — C 
diagram). However, a more natural choice, especially in systems with multiple companions, would 
be the system's barycenter. From Fig. [U if z measures the distance of the single binary object from 
the line perpendicular to the line of sight and passing through the systems barycenter, then the 
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eclipsing period variation (or the observed minus computed, O — C, timing difference) will be 



n 



c 



Kb, 



1 + eb^i cos fb,i 



sin(/b,i + LVb, 



(1) 



where Kb^i = ab^i sin Ib^i / c measures the semi-amplitude of the LITE orbit in the O — C diagram 
with c denoting the speed of light and i = 1, 2 referring to the LITE binary orbit when considering 
either one of the two companions. In equation (1), ab^i is the semi-major axis of the binary (single 
object) orbit, eb^i is its eccentricity, Ib^i is the orbital inclination with respect to the plane of the sky, 
ujb^i measures the binary's argument of pericenter and fb^i denotes its true anomaly. We note that 
we do not include any secular timi ng variation due to mass-transfer between the two components. 
As pointed out bv lLee et al.l (|201ll ). this effect is minimal and we here assume that it is negligible. 
Furthermore, it is important to note that i = 1,2 describes two separate two-body problems. For 
i = 1, we consider the binary and the inner companion whereas for i = 2, we consider the binary 
and the outer companion. To obtain the resulting period variation of the measured mid-eclipse 
times due to both companions, one then usually assumes the superposition principle and uses the 
expression r = ri + T2 for the combined light-travel time effect. 

The expression in Equation ([T]) is obtained by evaluating the z-component of the binary that 
is along the line of sight, in a co ordinate syst em with origin at the system's barycenter. This is 
different from the formulation in llrwini (j 19521 ) who employs a coordinate system with origin at 
the center of the LITE orbit. The result is the omission of the Cb^i sin Ub^i term wh en comparin g 
with the corresponding expression for the light-travel timing difference presented in llrwini (j 19521 ). 
We find that a barycentric coordinate system is more intuitive when carrying out the subsequent 
dynamical analysis. Since Tj is a quantity measuring a time difference, it follows that there should 
be no difference in the derived orbits when formulating the light-travel timing difference in the two 
co ordinate s y stems . We have tested the latter and used the best-fitting parameters from Table 6 
in Ilcc et al.l ( 2011 ) as an initial guess for the two formulations. The results are shown in Fig. [2] 
showing a shift along the secondary axis for the two best-fitting models. The final best-fitting 
or bital parameters were similar for the two models and were close to the parameters as determined 
bv lLee et all (|201lh . 



Finally, it is worth mentioning some properties of the LITE orbit for clarification. The previ- 
ously mentioned Kepler elements all describe the binary system's orbit (as a single object). None 
of these elements are directly attributed to the (possibly unseen) companions orbit. The Keple- 
rian elements of th e comp anion(s) are inferred only indirectly from first principles. We refer to 
Murray fc: DermottI (120011 ) for details of properties of orbits in a barycentric coordinate system. 
In the following we qualitatively outline some relationships between the two orbits and refer to 
Fig. [TJ First, the two semi-major axes are related to each other via the masses. Second, the orbital 
eccentricity, inclination and orbital periods (-Pi, 2) are the same for the single-binary object and 
the companions orbit. Also the apsidal lines of the two orbits are anti-aligned (see Fig. [1]) giving 
rise to a 180° difference between the two argument of pericenter angles. Orbital quantities related 
to a companion are denoted with a numeral subscript starting with the inner companion first (i.e 
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e^ shows t he ec centricity of the inner companion). We note that this convention is different in 



Lee et alj (|201ll ). In their work the authors use the subscript 4 (3) to denote the inner (outer) 



companion. 



3. STABILITY STUDY 



We examined the orbital stabihty of the proposed three-body system (binary and two M- 
type stars) of SZ Her. Orbital parameters of the two proposed c ompanions are listed in Table [T] 
with formal la uncertainties as reproduced from iLee et al.l (|201ll ). The masses and orbital semi- 
major axes are stated without formal uncertainties and were determined along with other orbital 
parameters, as outlined below. We show a graphical representation of the two osculating M-star 
orbits in Fig. [3] for two values of the argument of pericenter. The figures assume the case where 
Ib,i = h,2 = 90° with the combined binary pair placed at the origin, defining the dynamical 
center in an astrocentric system. An inspection of the derived osculating elements indicates that 
the pericenter distance ai^2(l — ^1.2) of the inner and outer companions are at 8.63 and 7.45 AU 
respectively, implying an orbit crossing geometry. Considering the large masses of the companions, 
such an orbital architecture is expected to be highly unstable. In the following, we will study the 
orbital stability of the two proposed circumbinary companions in more detail. 

We used the variable-step, Bu lirsch-Stoer (BS2) A-body a l gorithm as impl emented in the 
orbit integration package MERCUR'v[^ ( Chambers fc Migliorini 1997 : Chambers 1999). The code was 
compiled using the Intel Fortran Compilei0 on an Linux based platform equipped with an INTEL-i5 
(2.8 GHz) CPU. The initial time step was set to 0.01 days. During the integrations the maximum 
relative energy error was a few times 10~ 



-13 



In our calculations, we combined the masses of the two binary components and treated them 
as one single object. We used this simplifying assumption in order to b e consistent with the orbital 
parameters (and minimum masses) as derived from the LITE formalism ( Irwin|[l952 ) . We integrated 
the orbits of binary and companion bodies in a "binarycentric" system with the combined binary 
mass at rest (this system is also referred to as a non-inertial astrocentric frame). This requires a 
transformation of orbital elements from the coordinate system defining the LITE orbit (barycentric 
frame) to the binarycentric frame. Throughout this work, we consider the two companions to be 
on the same plane. That is, Ii = l2- This consideration may be justified noting that circumbinary 
objects may form in the same accretion disk where the binary system was formed, and as a result, 
orbits are more or less aligned with each other. In such a scenario, the inclination of the orbit of 
the system with respect to the plane of the sky is likely to be close to 90°. However, we remind the 
reader that no observational data exist for SZ Her that might allow the determination of Ii,2- 



www. arm. ac.uk/~jec 
*We used the following compiler flags at compile time: ifort -03 -real-size 64 
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3.1. INITIAL CONDITIONS FROM LITE ORBIT 



In general the inclination of the LITE orbit relative to the plane of the sky is not known. 
Only in the case of Ii^2 — 90° where the companion is also eclipsing one (or both) of the binary 
components information abou t this angle can be obtained f rom photorn e tric m easurements (e.g., 
Kepler-16b. loovle et all (j201ll ) and Kepler-34b, Kepler-35b. I Welsh et alJ (j2012l )). However, in the 
case of SZ Her, such a detection seems difficult due to the long LITE periods involved. Possi- 
bly third/fourth- light might be detectable from spectral data. The lack of constraints in orbital 
inclination implies that only information about the minimum mass (mi_2 sin/i_2) and minimum 
semi-major axis (ai^2 sinli^2) can be obtained. 

Information on the minimum mass from each individual LITE orbit is extracted somewhat sim- 
ilar to the case of a single-lined spectroscopic binary. The minimum mass of the two circumbinary 
companions were determined from a Newton-Raphson iteration using the individual massfunctions 

47r2(a6,jsinli)3 _ {rmsmlif 



GP? 



{Mhp + rrii 



(2) 



for each LITE orb it and the comb ined mass of the binary pair Mft„ ~ 2.32 Mq (|Lee et al.l 12011 



In agreement with lLee et al.l (j201ll ). the masses of the two M-dwarfs were found to be 0.19 Mq for 
the inner, and 0.22 Mq for the outer companion (see Table [T]). 

In calculating the values of the masses, we chose not to consider the formal uncertainties in 
the companions masses since as presented bv lLee et al.l (j201ll ). these (formal) errors are on a 0.1% 
level. From numerical experiments, these minute changes in mass have little impact on the overall 
dynamical evolution of the system. 



The minimum semi-major axes of the LITE orbits are adopted from lLee et al.l (j201ll ) and are 
stated in the barycentric frame. In order to determine the semi-major axes of the two companions 
in the astrocentric frame, we used Kepler's third law, since the minimum mass and orbital period 
along with the total mass of the binary system are known parameters. We found the minimum 
semi-major axis for the inner companion to be aisin/i = 16.6 AU. For the outer companion we 
found a2 sin/2 = 26.6 AU. 



3.2. RESULTS FROM ORBIT INTEGRATIONS 

In the following we present results from our stability study of the two proposed companions 
orbiting SZ Her. We carried out a series of short term integrations (up to 10^ years) for various 
initial conditions. We first considered the likely case of Ii^2 = 90° and assigned the companions 
their minimum masses, minimum semi- major axes and eccentricites. We then varied the initial 
argument of pericenter {001^2) and initial mean longitude (Ai_2) for the two companion M-type stars 
and integrated their orbits. Figure H] shows the results of a few integrations. The upper panel in 
this figure shows the case where the orbital apsidal lines of the two companions are initially parallel 
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and both companions start with the same mean longitude (A3^4 = 0°). The result is a quick capture 
of the two companions into a binary system. Although this initial condition results in a long-term 
stability, it does not reflect the derived Keplerian parameters of the two LITE orbits proposed by 



Lee et al.l (|201ll ) . The middle panel of Fig. [J] shows the result of an integration where the outer 
circumbinary component escapes the system. In this integration, we used the fitted argument of 
pericenter for each LITE orbit to derive the argument of pericenter for the two companions. Finally, 
we considered the case where the outer M-dwarf companion is started at its apocenter 02(1 — 62) 
with A2 = 180°. The remaining Kepler elements were chosen similar to the previous case. The 
corresponding orbits are shown in the lower panel of Fig. HI As shown here, for this set of initial 
conditions, the inner companion escapes from the system in less than 100 years. 

In general, for the escape scenarios, because of the conservation of energy, one companion is 
transfered to a smaller orbit as a result of the escape of the other companion. We encountered 
this behaviour in all our integrations where an escape took place. We then considered various 
orbits with semi-majo r axes and eccentricities within and beyond their formal la uncertainties 



as given by iLee et al.l ()201ll ). The majority of these integrations (especially those with higher 
orbital eccentricities) resulted in unstable orbits. Low-eccentricity orbits showed a somewhat longer 
lifetime. The main mechanism resulting in the system's instability was the escape of one of the 
companions as a result of a close encounter with the other. 

In order to search for a stable orbital configuration which reflects the main characteristic of 
the two proposed one-companion LITE orbits, we considered initial conditions with different LITE 
orbit inclinations for the two companions on co-planar orbits. Although the orbital separations 
increases with decreasing inclination, mutual gravitational perturbations are still important since 
the masses of the two companions increase as well. Therefore, it is not obvious apriori if a stable 
configuration can be found. We integrated the equations of motion for ten different values of the 
inclination with respect to the plane of the sky in the range of /i^2 £ [1°, 90°]. For each inclination, 
we scaled the minimum mass and semi-major axis for the two companions accordingly in order to 
obtain their true (assumed) mass. In each integration, we placed the outer component at A2 = 180°. 
A representative subsample of our results is shown in Fig. [5l In all ten cases we found the system 
to be highly unstable on short timescales. To double-check our results, we repeated the numerical 
integrations using the RADAU algorithm (also available in the MERCURY package). The outcome 
confirmed our previous results. 

In each integration, we found multiple close approaches that either led to ejections or collisions 
between the two proposed circumbinary companions. We repeated our integrations with various 
initial Ai^2 resulting in no difference in the overall inferred orbital instability. Considering a wide 
range of orbital possibilities, we find that all our orbital integrations to result in highly unstable 
systems. 
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METHODS AND RESULTS FROM REVISED LITE MODELS 



In order to find a possibly stable LITE model to the proposed quadruple system, we carried 
out an extensive search of the parar neter space. The a nalysis, methodology and technique are 
similar to those described in detail in iHinse et alj (j201ll ). The only exception is that we now 
fo rmulate the LITE model in the barycentric frame by omitting the e sin uj term in equation 3 
of iHinse et alJ (120111 ) as outlined earlier. In the following, we briefly repeat main aspects of the 
underlying analysis. 

We used the Leyenbe rg-Marquardt least-square minimisation algorithm as implemented in 
MPFIT (jMarkwardtl l2009l ) software. The goodness-of-fit statistic of each fit was evaluated from 
the weighted sum of squared errors, x^- Here we use the reduced chi-square statistic Xr- We seeded 
107,625 initial guesses within the framework of a Monte Carlo experiment. 

Each guess was allowed a maximum of 500 iterations before termination. Converged solutions 
were accepted with initial guess and final fitting parameters recorded. Initial gu esses of the mode l 



parameters were chosen at random from either a uniform or Gaussian distribution. iLee et al.l (j201ll ). 
for example, provided a Lomb-Scargle period analysis on the complete timing data set. They 
determined two possible dominant frequencies of 3.53 x 10"^ cycle and 7.06 x 10~^ cycle 
associated with the two proposed circumbinary companions. The shorter period is relatively well 
determined given the long observational time span of SZ Her. We, therefore, draw random periods 
from a Gaussian distribution with the standard deviation three times larger for the longer period. 
The standard deviation for the shorter period was set to ±8 years with its mean centered at 43 
years. This choice somewhat restricts the quasi- global exploration of the (in principle infinite) 
X^-parameter space to within a physically meaningful sub-region. 



We used the same timing data set as in 



Lee et al 



defined by a uniform time standard. In iLee et al.l (|201ll ). t iming data for SZ Her were recorded 



(j201ll ). but applied our model on times 



in the UTC time frame, which is known to be non- uniform (jCuinan &: Ribasll200ll ). We therefore 
transformed the HJD (Heliocentric Julian Date ) timing recor ds in UTC (Coordinated Universal 
Time) into the terrestrial time (TT) standard (jBastianl l2000l ) . The resulting unit of TT time is 
HJED (Helocentric Julian Ephemeris Date). 



4.1. Results 



m 



The result s frorn our Monte Carlo experiment are somewhat similar to the results presented 
Hinse et al.l (|201ll ). We discarded all guesses which reached a lower or upper boundary in 



one of the model parameters (since no formal errors are supplied for these parameters within the 
MPFIT algorithm). The total number of qualified guesses were then reduced to 30,700. The 
majority of guesses (15,659 or 50%) had final reduced goodness-of-fit paramet er in the i nterva l 



1.0086 < < 1.0186. This is in agreement with the two-LITE model found in iLee et al.l ( 201 
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However, a small number (3758 or 12%) of initial guesses had a final reduced goodness-of-fit statistic 
in the interval 0.9886 < Xr < 0.9986 with the best-fitting model resulting in Xr = 0.9886. All 
remaining fits had Xr > 1-05. 



Our best-fitting model is somewhat smaller than the two-LITE fit presented in iLee et al 



(j201ll ). We show the best-fitting model in Fig. [6] with the best fitting parameters listed in Table 
[2j We note that the orbital eccentricity o f the inner comp anion increased significantly from 0.48 
to 0.76 when comparing with the work in iLee et al.l ( 201ll ). Orbital radii, periods, as well as the 



minimum masses of the two companions are almost unchanged. We calculated the mean value and 
corresponding standard deviation of the final semi-major axis and eccentricity for the two LITE 
orbits. For the inner LITE orbit, the average final semi-major axis and eccentricity were 1.24zb0.38 
AU and 0.57 it 0.17, respectively. For the outer LITE orbit, the average final semi-major axis and 
eccentricity were 2.32 ± 0.41 AU and 0.59 ± 0.22. 



4.2. ORBIT STABILITY OF BEST-FIT AND LITE COMPUTATION 

In the previous section, we obtained an improved fit to the existing timing data of the eclipsing 
SZ Her system. The resulting best-fit osculating orbits of the two proposed companions showed an 
orbit-crossing architecture. To examine the orbital stability of the best-fit quadruple system the 
equations of motion were integrated in an astrocentric system using MERCURY. The binary pair was 
treated as a single object. We considered a large combination of different initial conditions and 
studied the final outcome for different values of semi-major axis, eccentricity and orbital angular 
variables. In particular, we studied the dynamics of the system near the suggested 2:1 mean- 
motion resonance. All integrations were carried out for 10000 years. Total system energy was 
conserved to within a few times 10"^^. In all cases, integrations resulted in the escape of one of 
the two companions from the system. We show four examples of the time evolution of the orbits 
of the two M-dwarfs in Fig. [71 It was assumed, initially, that the two companions were co-planar 
(/i^2 = 90°) with the binary system. The upper (lower) two panels of Fig. [7] show the case for 
wi = (cji = 180 degrees). In Fig. [7t, the outer companion escapes the system within a few years. 
The resulting system consists of a single companion on a stable Keplerian orbit causing (at most) 
a single sinusoidal LITE effect. 

To demonstrate this single sinusoidal variation, we computed the light-travel time effect from 
a direct numerical integration for two smaller values of the two companions masses. We considered 
the case for which the two companions are co-planar with the binary plane {Ii^2 = 90°). The 
upper panel of Fig. [8] shows an unstable system and is similar to the orbit shown in Fig. [7^. As 
explained earlier, one companion is ejected from the system within 50 years. The resulting LITE 
effect exhibits an initial variation in the binary period. Because of the conservation of the total 
linear momentum, the binary (single object) and bound companion move in a direction opposite 
to that of the ejected companion. The result is a one-component LITE effect superimposed on a 
constant period change (linear trend) due to the systemic motion of the whole system towards the 
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observer (negative O — C values). The systemic velocity is obtained from the slope of the linear 
trend. The lower panel of Fig. [8] shows the orbits of the two companions with their mass reduced by 
a factor of ten. The initial conditions in these simulations are similar to those in the upper panel of 
Fig. [8j Because of the smaller mass of the circumbinary companions, the system is now stable on 
a 200 year time scale. However, for longer integration times, the outer companion is ejected from 
the system as a result of strong mutual interactions with the inner one. 

Within the 200 year integration the resulting LITE effect shows a quasi-periodic modulation of 
the binary period. The semi-amplitude -fCfe^j of this signal is around 0.0015 days or 130 seconds. In 
order to reliably detect such a LITE signal would require la timing uncertainties to be significantly 
smaller than 112 seconds (0.0013 days). 



CONCLUSION & DISCUSSION 



In this work, we have reanalysed the complete timing data set of SZ Her and determined an 
improved two-LITE fit using an extensive Monte Carlo based Xr'P^'^^-^^ter space search. Our model 
used the center of mass of the one-companion system as the origin of the underlying coordinate 
system. Using this new approach, we find no significaii t chan ge in the derived parameters when 
comparing our results with those presented in lLee et al.l (j201ll ). As shown in this work, the LITE 
orbits derived from the two coordinate systems (one with origin in the center of the LITE orbit 
which includes esinw term and the other with origin in the system's barycenter as expressed by 
Equation ([1])) are similar and no significant difference was observed between the two models. 



In comparison with lLee et al.l ()201 



and assuming the two-LITE model is correct, our fit with 
= 0.989 provides a slightly better description of the underlying timing data. The existence of 
the improved fit presented in this work is most likely explained by our thorough (quasi-global) Xr 
parameter search. However, the majority of i nitial guess e s resu lted in a slightly larger Xr statistic. 



which is consist ent with the first fit found in iLee et al.l (j201ll ). To some extent, the fit found by 



Lee et al.l ()201ll ) is more "stability-friendly" due to the lower eccentricity of the inner proposed 
companion. Furthermore, our fit also suggests the two proposed companions to be in a near 2:1 
mean-motion resonance with periods of 42 and 90 years for the inner and outer orbits. However, 
we have considered various initial configurations close to this resonance all of which resulted in 
unstable orbits. 



Our stability analysis showed that all models presented in iLee et al.l (j201ll ) and in this work 
resulted in unstable quadruple systems (considering the three-body problem). In particular, our 
parameter study of the unknown orbital inclination with respect to the line of sight indicated that 
reducing the inclination from 90° to 1° results in an increase in the companions semi-major axis. 
However, at the same time, the companion's masses also increas which in turn introduces large 
gravitational perturbations between these two objects. Our inclination survey concluded that all 
low-inclination orbits also result in unstable systems on very short time scales. 
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If the observed timing variation is real and caused by the presence of two circumbinary objects, 
then the system needs to be in a stable configuration. We therefore conclude that either the two- 
LITE model is an inadequate description to the data, and/or the underlying data set is insufficient. 
For that reason, in order to constrain any future modeling efforts, we encourage the aquisition of 
additional photometric observations of SZ Her. 

There may also be an alternative, although unlikely, exp lanation for the obs erved period mod- 



ulation in the SZ Her binary orbit as originally outlined bv iHorner et al.l (j201ll ). This possibility 
depicts a scenario for which the two companions currently undergo a dramatic dynamical evolu- 
tion with a transition from a stable to an unstable configuration with one companion escaping the 
system shortly as suggested in this work. The dynamical reason for this instability would remain 
an open question and additionally renders this scenario unlikely. If true, the observational conse- 
quence of such a scenario should make it possible in the foreseeable future, to detect a linear trend 
in the measured O — C diagram. This was demonstrated numerically in the top panel of Fig. [8j 
However, due to the short orbital instability time scales found in this work, it seems unlikely that 
we are currently witnessing a break-up scenario in which the system enters the very endstate of its 
dynamical evolution with one companion on the verge of escaping the binary system. The leading 
question is: why should we be in the fortunate situation and observe the last few hundred years of a 
supposedly long-lived (possibly millions of years) system that gradually evolved towards a general 
instability? 

On the other hand, and as mentioned earlier, a dynamical inspection of currently known 
proposed circumbinary planets reveals these systems are unstable as well (at least three out of 
four). This finding could be attributed to the fact that ground-based photometric observations 
are more sensitive to detect sub-stellar circumbinary objects introducing large-amplitude period 
variations superimposed on the linear ephemeris of the mid-eclispe times of the binary. In this case 
the quintessence would be: larger masses would introduce larger perturbations and cause larger 
period modulations in the — C diagram and therefore such systems would be prone to disintegrate 
on a short time-scale due to strong mutual in teractions . Anot her argument for the non-existence of 



the two components with orbits proposed in lLee et al.l ( 20 111 ) comes from numerical iV-body orbit 



cal culations. Studies of the dynamical stability of hierarchical four-body systems were carried out 



by ISzell et al.l (120041 ). Their work suggests that low-mass stellar objects on circumbinary orbits 
result in unstable hierarchical systems. Considering symmetric pairs of masses, these authors 
showed t hat for large mass-r atios, the most likely event is a double binary configuration (e.g. 



HD98800, iFurlan et al.l (|2007l )). Indeed, in one of our experiments, we confirmed such an outcome 
from a direct numerical integration. On the other hand, for small mass-ratios (possibly comparable 
to planetary masses), their work suggest that the most likely outcome are circumbinary orbits. 
This again was demonstrated when decreasing the masses of the two companions by a factor of 
10. However, a thorough stability analysis of circumbinary four-body low-mass stellar systems, 
considering a large range of orbital parameters and masses, would be helpful to identify stable 
domains of circumbinary orbits. 
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It isimportant to note that other sources may exist that can create modulations in the ob- 
served binary period. This could be in the form of magnetic interaction between the two binary 
components, and/or mass or an gular moin e ntuin transfer resulting in a secular modulation of ob- 
served timings. In their study, iLee et al.l (j201lh point out the possibility that SZ Her (being a 
semi-detached binary pair with the less massive component filling its Roche lobe) is currently un- 
dergoing a phase of weak mass transfer. However these authors provide arguments that this effect 
is likely to be negligible. Other mechanism potentially capable of causing eclipse timing variations 
is the direct grav itational perturb ations by a companion on the binary orbit. This possibility was 
not considered in iLee et al.l (j201lh and is left for a future study. 



Star-spots also could introduce stellar jitter mimicing period variations as discussed in lWatson &: Dhillon 



(|2004l ). Also, timing measurements might suffer from systematic measurement err ors introducing 



correlated red noise possibly resulting in wrong model parameters (jColes et al.ll201ll ). Unaccounted 



systematic (red) timing errors resulted previously in the false detection of planets around pulsars. 

To unveil the true nature of the observed timing variation we encourage the aquisition of 
future photometric/spectroscopic follow-up observations of SZ Her allowing to further constrain 
and refine timing models. Future models favouring a two-companion solution should be tested for 
orbital stability and the resulting O — C variation obtained from numerical integrations compared 
with the inferred timing measurements. 
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Parameter SZ Her(AB)C SZ Her(AB)D Unit 

minimum semi-major axis, ai^2 sin /i^2 16.6 26.6 AU 

eccentricity, ei,2 0.48 ±0.17 0.72 ±0.09 

argument of pericenter, uji^2 105 ± 10 268.6 ± 7.5 degrees 

orbital period, Pi,2 42.5 ± 1.1 85.8 ±1.0 year 

minimum mass, mi^2 sin /i_2 0.188 0.222 Mq 

Table 1: Binarycentric or astrometric orbital parameters of the two circumbinary companions C 
and D. Note that we have accounted for the 180° d ifference in the a rgument of pericenter angle. 



Parameters with formal la uncertainties are from iLee et al.l (j201ll ). All other parameters are 
obtained as outlined in the text. 



- 16 - 



Parameter 


two-LITE 


Unit 




n (i = 1) 


r2 {i = 2) 






0.989 


- 


RMS 


0.00339 


days 


To 


2,434,987.38455(31) 


HJED 


Po 


0.818095801(11) 


days 


ab^i sin Ib,i 


1.49 ±0.53 


2.23 ±0.41 


AU 


eb,i (or 61,2) 


0.76 ±0.14 


0.69 ±0.12 






2.42 ±7.1 


286.16 ± 10.1 


deg 


Tb^i 


2,422,649(322) 


2,439,551(198) 


HJED 


Pb,^ (orPi,2) 


15286 ± 224 


32860 ± 341 


days 


Kb,i 


0.00864(32) 


0.00129(22) 


days 


fim) 


0.00191(65) 


0.00138(37) 




rrii sin /j 


0.23 ±0.08 


0.21 ±0.05 


Mq 


Qi sin li 


16.5(9) 


27.4(8) 


AU 


&% 


0.76 ±0.14 


0.69 ±0.12 






182 ±7.1 


106 ± 10.1 


deg 


P^ 


15286 ± 224 


32860 ± 341 


days 



Table 2: Best-fit parameters for the LITE orbits of SZ Her corresponding to Fig. [6l Subscripts 1, 2 
refer to the circumbinary companions with i = 1, the inner, and i = 2, the outer, companions. 
RMS measures the root-mean-square scatter of the data around the best fit. Formal uncertainties 
obtained from the covariance matrix are valid for the last digit and shown in paranthesis. Note 
that the eccentricity and orbital period are shared quantities as outlined in the text. The last five 
lines are quantities of the two companions in the astrocentric coordinate system. 
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Fig. 1. — Graphical outline of the LITE orbit for the case Ih^i = 90°. The point O is the binary- 
companion center of mass and denotes the origin of the coordinate system. The x-axis points 
towards the line of reference (T). The (x,?/)-plane coincides with the plane of the sky and is 
perpendicular to the line of sight. The LITE binary orbit is shown as a solid ellipse with the 
instantaneous position of the binary at P. The angles to and / denote the argument of pericenter 
and true anomaly, respectively. 
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Fig. 2. — Comparing best-fit LITE models in the ellipsoid centered (upper panel with Xr — 1-0129) 
and barycentric centered (bottom panel with = 1 .0127) coordinate systems. The initial guesses 
resulting in the two fits were taken from iLee et al.l (j201ll ). In each panel, the inner (outer) com- 
panion is denoted by ti (r2). In both panels the root-mean-square (RMS) scatter of data points 
around the best fit is around 300 seconds. 
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Fig. 3. — Geometry of the two orbits (in the orbital plane) corresponding to the LITE fit parameters 
in Tabled! Both companion orbits are relative to the binarycentric reference frame (with the binary 
taken to be a single object) with origin at the center of the cross-hair. The orbits were integrated 
for one orbital period considering massless objects to visualise the two osculating single LITE 
orbits. We show both the w = 0° (la, 2a) orbits and the orbits for u = (285 — 180)° (lb) and 
oj = (88.6 - 180)° (2b). 
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Fig. 4. — Demonstrating unstable orbital time evolution (in the orbital plane) of the two companions 
for different initial conditions. The stipulated orbit represents the outer companion. The center of 
the cross-hair marks the origin of the binarycentric reference frame. In all cases we consider /i_2 = 
90°. Upper panel: a;i,2 = 0°,Ai,2 = 0°. Middle panel: ui = (286-180)°,a;2 = (89-180)°, Ai,2 = 0°. 
Bottom panel: wi = (286 - 180)°, ^2 = (89 - 180)°, Ai = 0°,A2 = 180° 
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Fig. 5. — Unstable orbits from numerical integrations for various LITE orbital inclinations. Each 
panel shows the orbits (in the orbital plane) using initial conditions that considers a scaled semi- 
major axis and mass for the two companions simultaneously. The inclination between the line of 
sight and the plane of the sky were Ii,2 = 5°, 10°, 30°, 50°, 70°, 80°. In ah integrations, A2 = 180°. 
The orbits of the two companion were assumed to be co-planar. The stipulated line always repre- 
sents the outer binary companion. The centre of the hair-cross marks the origin of the barycentric 
reference frame. 
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Fig. 6. — Best fit model with Xr = 0.989 (dash-dot-dot-dot) from our many-guess Monte Carlo 
experiment and the two companion sinusoidal-like variations: inner (dash) and outer companion 
(dash-dot). The corresponding orbital parameters for the two LITE orbits are shown in Table [2j 
The root-mean-square scatter of data around the best fit is ~ 292 seconds. We show the three 
zblo" timing error bar s in th e lower left corner corresponding to 0.0036, 0.0020 and 0.0013 days as 
adopted by iLee et al.l (j201ll ) for various observation techniques. 
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Fig. 7. — Time evolution of the two M-dwarf companions using initial conditions from Table EJ 
The mean longitude is denoted by A. Upper two panels are for wi = and bottom two panels for 
ui = 180°. In each panel the binary pair is placed at the origin of the coordinate system. 
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Fig. 8. — Numerical computation of the orbit (left) and the resulting LITE effect (right) for two 
different scenarios of companion mass. Initial conditions are the same as in Fig. [7^ with Ii^2 = 90°. 
Upper panel: Companions masses are mi = 0.23 Mq and m2 = 0.21 Mq. Lower panel: Companions 
with masses of mi = 0.023 Mq and m2 = 0.021 Mq. Vertical bars in the right panels represent 
±1(7 uncertainties with a corresponding to 0.0013 (112 seconds), 0.0020 (173 seconds) and 0.0036 
days (311 seconds). 



